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^ ; Abstract 

The spectrum of cold dark matter particles on Earth is expected to have 
> ■ peaks in velocity space associated with particles which are falling onto the 

Galaxy for the first time and with particles which have fallen in and out of 
the Galaxy only a small number of times in the past. We obtain estimates 
for the velocity magnitudes and the local densities of the particles in these 
5_i ■ peaks. To this end we use the secondary infall model of galactic halo formation 

which we have generalized to take account of the angular momentum of the 
dark matter particles. The new model is still spherically symmetric and it 
admits self-similar solutions. In the absence of angular momentum, the model 
produces flat rotation curves for a large range of values of a parameter e which 
is related to the spectrum of primordial density perturbations. We find that 
the presence of angular momentum produces an effective core radius, i.e. it 
makes the contribution of the halo to the rotation curve go to zero at zero 
radius. The model provides a detailed description of the large scale properties 
of galactic halos including their density profiles, their extent and total mass. 
We obtain predictions for the kinetic energies of the particles in the velocity 
peaks and estimates for their local densities as functions of the amount of 
angular momentum, the age of the universe and e. 
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I. INTRODUCTION 



Experiments are presently under way which attempt to identify the nature of dark matter 
|]J by direct detection on Earth. The dark matter candidates which are being searched for in 
this manner are WIMPs and axions. WIMPs is an acronym for "weakly interactive massive 
particles" . The best motivated candidate of this type is the lightest supersymmetric partner 
in supersymmetric extensions of the standard model of particle physics || . The mass range 
for which WIMPs provide the critical energy density for closing the universe is a few GeV 
to a few hundred GeV. The axion is a light pseudo-scalar particle whose existence has been 
postulated to explain why, in the context of the standard model of particle physics, the 
strong interactions conserve P and CP [[|. The likely mass range for which axions provide 
the critical energy density for closing the universe is 10~ 4 eV < m a < 10~ 7 eV Q. 

Axions and WIMPs are the leading cold dark matter (CDM) candidates. Other forms 
of dark matter are neutrinos and dark baryons. From the point of view of galaxy formation, 
the defining properties of CDM are: 

1. that CDM particles, unlike baryons, are guaranteed to interact with their surroundings 
only through gravity, and 

2. that CDM particles, unlike neutrinos, have negligibly small primordial velocity disper- 
sion. 

Studies of large scale structure formation support the view that the dominant fraction of 
dark matter is CDM. Moreover, if some fraction of the dark matter is CDM, it necessarily 
contributes to galactic halos by falling into the gravitational wells of galaxies and hence 
is susceptible to direct detection on Earth. WIMPs are being searched for by looking for 
WIMP + nucleus elastic scattering in a laboratory detector ||. The nuclear recoil can be put 
into evidence by low temperature calorimetry, by ionization detection or by the detection of 
ballistic phonons. Axions are being searched for by stimulating their conversion to photons 
in a laboratory magnetic field f||7]]. The experimental apparatus involves an electromagnetic 
cavity placed in the bore of a superconducting solenoid. When the resonant frequency of 
the lowest TM mode of the cavity equals the axion mass ihv = m a c 2 ), some galactic halo 
axions convert to microwave photons inside the cavity. If a signal is found in the cavity 
detector of dark matter axions, it will be possible to measure the energy spectrum of the 
axions with great precision and resolution. (The energy resolution is limited only by the 
measurement integration time and the stability of the local oscillator with which the axion 
signal is compared.) Hence the question arises: what can be learned about our galaxy and 
the universe by analyzing such a signal? The main purpose of this paper is to address this 
question by predicting properties of the CDM spectrum on Earth in terms of cosmological 
input parameters. Incidentally, all CDM candidates have the same phase space distribution, 
and hence the same spectrum on Earth, in the limit where their small primordial velocity 
dispersions are neglected. We are motivated in part by the fact that knowledge of the 
spectrum may help in the discovery of a signal. 

In many past discussions of dark matter detection on Earth, it has been assumed that 
the dark matter particles have an isothermal distribution, or an adiabatic deformation of 
an isothermal distribution. A strong argument in favor of this assumption is the fact that 
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it predicts the rotation curve of the galaxy to be flat at large radial distances. Indeed, self- 
gravitating isothermal spheres always have density distributions p(r) which fall off at large 
r as 1/r 2 . Moreover, they have a "core radius", i.e. a radius a within which the density p(r) 
no longer behaves as 1/r 2 but goes to a constant as r — ► 0. The behaviour may, for most 
practical purposes, be approximated by the function p(r) = p(0)[l + (r/a) 2 ] -1 . Thus, the 
contribution of an isothermal halo to the galactic rotation velocity goes to zero as r — > 0. 
This feature of isothermal halos is attractive as well because it is known that, in spiral 
galaxies like our own, the rotation velocity at small radii may be entirely accounted for by 
the bulge and the disk. In our galaxy, the core radius is such that roughly half of the rotation 
velocity squared at the solar radius r s ~ 8.5 kpc is due to the disk and bulge while the other 
half is due to the dark halo. Thermalization of the galactic halo has been argued to be the 
outcome of a period of "violent relaxation" || following the collapse of the protogalaxy. If 
it is strictly true that the velocity distribution of the dark matter particles is isothermal, 
then the only information that can be gained from its observation is the corresponding virial 
velocity and our own velocity relative to its standard of rest. 

However, one may convince oneself that the velocity distribution of dark matter particles 
has a non-thermal component. Consider the fact that our closest neighbor on the galactic 
scale, the galaxy M31 in Andromeda, at a distance of order 730 kpc from us, is falling towards 
our galaxy with a line-of-sight velocity of 120 km/sec. This motion can be understood to be 
due to the mutual gravitational attraction between the two galaxies: first they were receding 
from each other as part of the general Hubble flow, this relative motion was halted and now 
they are falling towards one another. We may use M31 as an indicator of the motion of any 
matter in our neighborhood. Moreover, if cold dark matter exists, then there is cold dark 
matter at every physical point in space (including everywhere we see nothing and which 
appears empty), because by Liouville's theorem the 3-dim. sheet in 6-dim. phase-space on 
which the CDM particles lie can not be ruptured. The thickness of that sheet is the tiny 
primordial velocity dispersion of the CDM particles, of order 10 -12 for WIMPs and 10~ 17 for 
axions (c = 1). The implication of the above is that, if CDM exists, there are CDM particles 
falling onto our galaxy continuously and from all directions. The motion of these particles 
gets randomized by gravitational scattering off giant molecular clouds, globular clusters and 
other inhomogeneities but complete thermalization of their velocity distribution occurs only 
after they have fallen in and out of the galaxy many times. As a result there are peaks in the 
velocity distribution of CDM particles at any physical point in the galaxy B. One peak is 
due to particles falling onto the Galaxy for the first time, one peak is due to particles falling 
out of the Galaxy for the first time, one peak is due to particles falling in for the second 
time, and so on. In particular, this is true on Earth. The width of the first two peaks, which 
we label n=l, is related to and is of order the velocity dispersion of the particles before they 
fall in for the first time. The width of the next two peaks (n=2) is expected to be somewhat 
larger as a result of scattering of the particles off inhomogeneities in the galaxy. The width 
of the next two peaks (n=3) is larger still because these particles have been scattered more. 
And so on. 

One of the main purposes of this paper is to obtain estimates of the sizes and velocity 
magnitudes of the velocity peaks on Earth. By "size" , we mean the contribution of the peak 
to the local mass density of the halo. By "velocity magnitude" , we mean the magnitude of 
the velocity vector of the particles in the peak as measured in the rest frame of the Galaxy, 
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i.e. in a frame which is not corotating with the disk. The tool we use is the secondary 



infall model [10] of galactic halo formation. This model assumes a single overdensity in an 



expanding universe. A halo forms around the overdensity because dark matter keeps falling 
onto it. The dark matter is assumed to have gravitational interactions only and to have 
zero initial velocity dispersion. The model also assumes spherical symmetry. Finally, in 
its original form, it assumes that the dark matter particles have zero angular momentum 
with respect to the center and hence that their motion is purely radial. Much progress 
Tl"l , |T2"ll in the analysis of the model was made as a result of the realization that the time- 



evolution is self-similar provided Q = 1 and provided the initial overdensity has a special 
scale-free form; see Eq.( 3.11|) . The parameter e that appears in this ansatz is related to 
the slope of the power spectrum of primordial density perturbations at the galactic scale. 
Self-similarity means that the halo is time-independent after all distances are rescaled by 
an overall time-dependent size R(t) and all masses by a time-dependent mass M(t). The 
self-similar solutions can be obtained numerically with great precision and some of their 
properties may be derived analytically. When the parameter e is in the range < e < 2/3, 
the density profile p(r) ~ 1/r 2 and thus the rotation curve is flat. 

However, for the purpose of estimating the sizes of velocity peaks, the secondary infall 
model without angular momentum is rather inadequate. In particular, it tends to overesti- 
mate the size of the peaks due to dark matter particles falling in and out of the galaxy for 
the first time. Indeed, angular momentum has the effect of keeping infalling dark matter 
away from the galactic center and this effect is largest for particles falling into the galaxy 
last. On the other hand, the presence of angular momentum destroys spherical symmetry 
and thus makes the actual evolution far more complicated and untractable. However, as 
will be explained in detail below, it is possible to include the effect of angular momentum 
into the secondary infall model without destroying its spherical symmetry by averaging over 
all possible orientations of an actual physical halo ||13|| . Moreover, the time evolution of 
the model with angular momentum thus included is still self-similar provided the angular 
momentum distribution is of a particular scale-free form. It was also found JT3] that angular 
momentum has the effect of making the halo contribution to the galactic rotation curve go 
to zero at the galactic center, thus introducing an effective core radius for the halo mass 
distribution. We define the effective core radius b to be the radius at which the halo con- 
tributes half of the galaxy's rotation velocity squared. For our Galaxy, b is of order our own 
distance to the Galactic center. This by itself suggests that the effect of angular momentum 
on the velocity peaks on Earth is not small. The model with angular momentum can be 
accurately solved on a computer. Its predictions for the effective core radius b, the local 
halo density, and the expected sizes and velocity magnitudes of the first few velocity peaks 
are tabulated below for representative values of the input parameters, which are the age 
of the Universe, the parameter e and the average amount of angular momentum. We also 
give an analytical treatment of the model under simplifying but realistic assumptions. It 
yields general formulae which may be used to estimate the expected sizes and the velocity 
magnitudes of the velocity peaks for a wide range of the input parameter values. 

In Section II we review the arguments of ref |9j why velocity peaks in the cold dark matter 
spectrum on Earth are expected, and add some comments of our own. In Section III we give 
a detailed description of the self-similar infall model, without and with angular momentum. 
In Section IV we describe how some of the model parameters are determined in terms of 
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observed properties of our Galaxy and we give the results of the numerical integration of the 
model. Section V contains our analytical treatment of the model. Section VI summarizes 
our results. 



II. PHASE SPACE STRUCTURE OF COLD DARK MATTER HALOS. 

In cold dark matter scenarios, the initial phase space distribution is a very thin sheet 
near v = Hf, where H is the Hubble rate and r is the position relative to an arbitrarily 
chosen reference point. The deviations from perfect Hubble flow which are present are 
associated with the primordial density perturbations that will produce galaxies and large 
scale structures by gravitational instability. Where a galaxy (or some other object) forms 
and grows, the phase space sheet is folding itself up. The process is illustrated in Fig. [I] 
for the simplified case where a single spherically symmetric overdensity is present in an 
otherwise homogeneous universe and where all dark matter particles move on radial orbits 
through the center of the overdensity. The line in the figure indicates the location of the 
dark matter particles in (r, f) phase space at an instant of time, r is the distance to center 
of the overdensity and f = dr/dt. As time goes on the line "winds up" in the clockwise 
direction, rotating most rapidly at the center. 

Fig. |l] shows that the velocity spectrum of cold dark matter particles on Earth, or 
anywhere else in the galaxy, has a series of peaks. One peak is due to particles falling onto 
the galaxy for the first time, passing by Earth while going towards the galactic center. A 
second peak is due to particles which are falling out of the galaxy for the first time, passing 
by Earth while going away from the galactic center. A third peak is due to particles falling 
onto the galaxy for the second time. A fourth peak is due to particles falling out of the 
galaxy for the second time. And so on. A rough estimate of the number N of velocity peaks 
on Earth in this idealized case may be obtained as the ratio of the age of the galaxy (~ 10 10 
years) to the time (~ 0.5 x 10 8 years) it takes a particle to fall to the center of the galaxy 
starting from rest at the Earth's location, with the result iV ~ 200. However, the presence 
of angular momentum of the dark matter particles tends to decrease N by restricting the 
range of radii over which dark matter orbits vary. (In the extreme limit of circular orbits, 
N = 1.) As will be seen below, this expectation |J is confirmed by our calculations. We 
shall also find that the number of peaks depends upon e. 

Of course the description of a galactic halo presented in Fig. [1] is much simplified. In the 
remainder of this section, we discuss the sensitivity of the conclusion, that there are peaks 
in the cold dark matter velocity distribution on Earth, to the simplifying assumptions that 
were made. In particular, we inquire into the effect of 

1. the gravitational scattering of the dark matter particles by inhomogeneities in the 
galaxy 

2. the angular momentum that the dark matter particles have with respect to the galactic 
center 

3. the velocity dispersion that the dark matter particles have before they fall onto the 
galaxy. 
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A. Scattering by inhomogeneities in the galaxy. 



The effect of the gravitational scattering of the dark matter particles by the inhomo- 
geneities in the galaxy, such as stars, globular clusters and large molecular clouds, is to 
"fuzz up" the phase space sheets. Consider a phase space sheet which in the absence of 
scattering produces on Earth a stream with unique velocity v. The collective effect of scat- 
tering by a class of objects of mass M and density n is to diffuse the velocities in the stream 
over a cone of opening angle AO given by 



where the time integral is over the past history of the particles in the stream, b is the impact 
parameter of a scattering and v is the velocity of the sheet relative to the scattering center. 
In the galactic disk, the giant molecular clouds are most likely the main contributors. With 
M ~ 1O 6 M , n ~ 3 kpc -3 , b max ~ kpc and b min ~ 20 pc, they yield AO ~ 0.05 for dark 
matter particles that have spent most of their past in the galactic disk. The contributions 
due to globular clusters (M ~ 5x 10 5 M Q , n ~ 0.3 kpc" 3 ) and stars (M ~ M , n ~ 0.1 pc" 3 ) 
are less important. At any rate, peaks due to dark matter particles that have spent much 
of their past in the central parts of the Galaxy are likely to be washed out. On the other 
hand, the peaks due to dark matter particles which have fallen in and out of the Galaxy 
only a small number of times in the past are not erased by scattering. 



The presence of the rotating galactic disk clearly indicates that the baryons in our galaxy 
carry angular momentum. This angular momentum is thought to have been produced by 
the gravitational forces of nearby galaxies when ours started to form. One should expect the 
dark matter in the galactic halo to have similar amounts of angular momentum and hence 
to move on non-radial orbits. If an infalling particle's angular momentum is large enough, 
its distance of closest approach to the galactic center is larger than our own distance (~ 8.5 
kpc) to the galactic center and hence it can not possibly reach us. It is nonetheless true that 
particles falling onto the galaxy for the first time reach us at all times, even if the typical 
distance of closest approach of such particles to the galactic center is much larger than 8.5 
kpc. 

Indeed, consider all particles that reach their turnaround radius at a given time. 
"Turnaround" refers to the moment in a particle's history when it reaches zero radial veloc- 
ity with respect to the galactic center for the first time, after its initial Hubble flow velocity 
has been halted by the gravitational pull of the galaxy and before it starts to fall onto the 
galaxy for the first time; see Fig. [l|. All particles that reach their turnaround radius at a 
given time are on a surface which, from a topological viewpoint, is a sphere enclosing the 
galactic center. Let us call this surface the "turnaround sphere" . By consulting a catalog of 




n (300km/s) 3 
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B. Angular momentum. 
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the galaxies in our neighborhood and plotting their radial velocities as a function of distance, 
one concludes that the radius of the present turnaround sphere is of order 1-2 Mpc for our 
galaxy. Consider the turnaround sphere at an arbitrary time t. At any point on that surface 
the angular momentum vector has a unique value parallel to the surface. Now, it is well 
known that a continuous vector field on a 2-sphere can not everywhere differ from zero. It 
must have at least two zeros. 

Hence there are two places on any turnaround sphere where the angular momentum 
vanishes. The particles from these two locations will pass through the galactic center when 
they fall onto the galaxy next, producing two velocity peaks there. By continuity (the 
phase space sheet can not tear), other particles on the turn-around sphere will produce two 
velocity peaks at any point sufficiently close to the center. Fig. |2| shows the time evolution of 
a turn-around sphere which is initially rigidly rotating about an axis and which subsequently 
is moving under the influence of the gravitational potential of an isothermal sphere. The 
figure demonstrates that any turnaround sphere passes (at least) twice by any point inside 
of it, once on the way in and once on the way out, assuming only that the point is inside 
the sphere both at its first and its second turnaround. By definition, second turnaround is 
when the sphere reaches its maximum size for the second time in its history, just after its 
first oscillation. Thus we find that when angular momentum is included, there are still two 
(possibly more but necessarily an even number) velocity peaks on Earth due to particles 
falling through the galaxy for the first time, two peaks or more due to particles falling 
through the galaxy for the second time, and so on. As we saw in the preceding subsection, 
these peaks are not erased by scattering off stars, globular clusters or giant molecular clouds. 
The sizes and velocity magnitudes of these peaks constitute the main topic of this paper. 

Finally, let us note the effect angular momentum has on the caustics of a halo. A 
"caustic" is a place where the dark matter density is large because the phase space sheet folds 
back there. The density actually diverges at the caustic in the limit where the thickness of 
the phase space sheet goes to zero. There is an outer caustic surface near the n-th turnaround 
radius with n = 2,3,4... ; see Fig. [H It can be shown that near a caustic surface, the dark 
matter density behaves as p ~ Q(x)/yfx + constant, where x is the distance to the caustic 
surface and Q(x) is the Heaviside function: 0(x) = 1 for x > and Q(x) = otherwise. 
Now, when angular momentum is absent, the center of the galaxy is a very special point 
because all dark matter particles go through the center at each orbital oscillation. The dark 
matter density goes as p ~ 1/r 2 at the center if there is no angular momentum and the 
rotation curve is flat. Thus, in the absence of angular momentum the galactic center is a 
caustic point. When angular momentum is present, that caustic point spreads into a set of 
inner caustic rings. Fig.|2] shows the appearance of such an inner caustic ring for the case of 
axial symmetry. The fact that the caustic appears has nothing to do with axial symmetry 
however. Rather, it is a consequence of the fact that when a sphere is turned "inside out", 
as illustrated in Fig.|^, a ring singularity must appear on the surface at some point during 
the process. Generally, the caustic ring appears near the place where the particles with 
the most angular momentum on a given turn-around sphere turn back at their distance of 
closest approach to the galactic center. Outer caustics at the Earth's location are likely 
to be very much degraded by scattering of the dark matter particles off inhomogeneities in 
the galaxy. However, inner caustics associated with particles which have gone through the 
central parts of the galaxy only a small number of times in the past are not much degraded. 
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The dark matter density on Earth could be much enhanced if we happen to be close to an 
inner caustic. 



C. On the velocity dispersion of infalling cold dark matter. 

In this study, when obtaining estimates of the average sizes and of the velocity magnitudes 
of the velocity peaks, we neglect the velocity dispersion 5vi n the cold dark matter has when 
it falls onto the galaxy for the first time. Presumably, this is a valid approximation provided 
5vi n is much smaller than the velocity dispersion 5v ga i ~ 10~ 3 of the galaxy as a whole. We 
argue in this section that this condition is probably satisfied although we will not attempt 
to provide a reliable estimate for the size of Sv in . The width 5v n of the velocity peaks due to 
particles falling in and out of the galaxy for the n-th time, where n is sufficiently small that 
the broadening effect of scattering of the particles by the galaxy's inhomogeneities can be 
neglected, is related to 5v in by Liouville's theorem: 5v n = 5fi n (t* jn )[p n /p(t ;( , in )] 1 / 3 where p n 
is the contribution to the local halo density from particles in the n-th peak, and 5vi n (t* tn ) 
and p(t*, n ) are the velocity dispersion and density those particles had at the time £* >n of 
their first turn-around. 

Let us emphasize that the values of the peak widths 8v n may some day be measured in a 
direct CDM detection experiment and that such data would provide information about our 
universe which is not readily accessible by other means. Also, if the widths of some peaks 
are small enough the sensitivity of the cavity detector of dark matter axions is improved 
by looking for narrow peaks. In the case of the present LLNL experiment ]7), which does 
look for narrow peaks in addition to looking for a signal whose width is set by the galaxy's 
overall velocity dispersion 5v ga i ~ 10 -3 , the sensitivity of the search is improved if there is 
a velocity peak with 5v n less than about 10 -8 and with a fraction of the local density larger 
than about 1%. 

Turning to the question how large 5vi n may be, let us start by describing the primor- 
dial velocity dispersion which is the contribution that is present even if the universe were 
completely homogeneous, i.e., it is the value of 5v in if our galaxy were the only density 
perturbation in the universe. For WIMPs, the primordial velocity dispersiom 5vw is due 
to the finite temperature Tb the WIMPs have when their kinetic energies decouple from 
the primordial heat bath. Thus 5vw ~ (2Tb / 'm) 1 ^ 2 (Rd / Rq) , where m is the WIMP mass, 
and Rd and Rq are the scale factors at temperature Tb and now. For m ~ 50 GeV and 
Tb ~ 1 MeV, one has v w ~ 10 -12 which is very small. For axions, the primordial velocity 
dispersion is due to the inhomogeneity of the axion field at temperature T\ ~ 1 GeV and 
time t\ ~ 2 x sec when the axion mass becomes equal to the Hubble expansion rate. If 
there is no inflation after the Peccei-Quinn phase transition at which the Upq(1) symmetry 
gets spontaneously broken, then the scale of inhomogeneity of the axion field is of order 
the horizon scale at time t\ and the primordial axion velocity dispersion today is therefore 
Sv a ~ (m a ti) _1 (_Ri/i? ) ~ 10~ 17 x (10~ 5 eV/m a ). If there is inflation after the Peccei-Quinn 
phase transition, then the axion field gets homogenized over enormous distances and 5v a is 
exponentially small. 

The primordial velocity dispersion, 5vw or 5v a , discussed in the preceding paragraph is 
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the thickness of the CDM phase space sheet. It constitutes a lower bound on the velocity 
dispersion 5vi n of infalling CDM. Additional velocity dispersion is expected because the 
phase space sheet may wrap itself up on smaller scales than that of the galaxy as a whole, 
as illustrated in Fig|| The phase space sheet wraps itself up wherever an overdensity has 
grown by gravitational instability past the linear regime {Sp/p < 1) into the non-linear one 
{Sp/p > 1). In theories of structure formation based upon cold dark matter, the spectrum of 
primordial density perturbations is flat, i.e., it has approximately equal power on all length 
scales. The matter density perturbations do not grow till the time t eq of equality between 
the matter and radiation energy densities. After t eq , all density perturbations which have 
wavelength less than the horizon (this includes all length scales of order a galaxy size) grow 
together at the same rate and therefore they all reach the non-linear regime at approximately 
the same time. In the standard CDM cosmology, the smaller scale clumps reach the non- 
linear regime somewhat earlier because the processed spectrum of density perturbations is 
not exactly flat on galaxy scales but is slightly tilted with more power on small scales. 

What happens in the non-linear regime is far from obvious. The rate of growth of an 
overdensity in the non-linear regime is of order ^Gp, where p is its mean density. Indeed 
1/y/Gp is of order the free infall time, which is also the time necessary to produce a new 
fold in the phase space sheet. At the start of the non-linear regime, as we just argued, 
all overdensities have densities of the same order of magnitude and they therefore grow at 
comparable rates by locally wrapping up the phase space sheet. However, overdensities of 
large physical size will tidally disrupt and therfore inhibit the growth of overdensities of 
smaller physical size in their neighborhood. In the immediate vicinity of our galaxy, there 
are no visible overdensities other than nearby dwarf galaxies such as the Magellanic clouds. 
Dark matter and accompanying baryons are nonetheless falling onto the galaxy now for the 
first time. It is the velocity dispersion of this unseen matter that we are interested in. If this 
matter is in large clumps, one might expect it to light up stars and thus become visible. On 
the other hand, it could be in clumps which have not lit up for some reason. However, any 
known object smaller than a galaxy (e.g. stars, globular clusters, large molecular clouds, 
dwarf galaxies) has velocity dispersion smaller than 5v ga i ~ 10~ 3 , and dark matter objects 
should be expected to be less clumped than baryonic objects because they can not dissipate 
their energy. On this basis, it seems safe to assume that 5vi n is considerably less than Sv ga i. 

There is a particular kind of dumpiness which is expected to affect axion dark matter if 
there is no inflation after the Peccei-Quinn phase transition. This is due to the fact that cold 
dark matter axions are inhomogeneous with Sp/p ~ 1 over the horizon scale at temperature 
Ti ~ 1 GeV when they are produced at the start of the QCD phase-transition, combined 
with the fact that their velocities are so small that they do not erase these inhomogeneities 
by free-streaming before the time t eq when matter perturbations can start to grow. These 
particular inhomogeneities in the axion dark matter are immediately in the non-linear regime 
after time t eq and thus form clumps, called "axion mini-clusters" flll|-|T6||. These have (16| 
mass M mc ~ 10 _13 M Q and size l mc ~ 10 12 cm, and therefore their associated velocity 
dispersion v mc = ^GM mc /l mc ~ 1CT 10 at time t eq . This velocity dispersion increases by 
about a factor 10 from t eq till the onset of galaxy formation because of the hierarchical 
clustering of the axion mini-clusters. This yields v mc ~ 10~ 9 as the contribution of mini- 
clusters to the velocity dispersion 5v j n of infalling axions if there is no inflation after the 
Peccei-Quinn transition. 
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III. SELF-SIMILAR INFALL MODELS. 



A. The radial infall model. 



The tool we use to obtain estimates of the sizes and the velocity magnitudes of the 
highest energy peaks is the secondary infall model of galactic halo formation. In its original 
form, this model is based on the following assumptions: 

1. the dark matter is non-dissipative 

2. it has negligible initial velocity dispersion 

3. the gravitational potential of the galaxy is spherically symmetric and is dominated by 
the dark matter contribution 

4. the dark matter particles move on radial orbits through the galactic center. 

Assumption 1 means that the only force acting upon the dark matter particles is the grav- 
itational pull of the galaxy. Assumption 2 states that before the galaxy starts to form, at 
some initial time tj, all the dark matter particles at the same position r*j relative to the 
galactic center move with the same velocity Vi. Provided tj is chosen early enough, this 
initial velocity is given by the Hubble expansion: 



H(ti) is the Hubble rate at time tj. The issue of the validity of assumption 2 is discussed at 
length in the previous section. Henceforth, we will take its validity for granted. Assumption 
3 is realistic because the gravitational potential of a galaxy as a whole (luminous plus dark 
matter) is nearly spherically symmetric even if the distribution of its luminous matter is not 
spherically symmetric at all. Assumption 4 is the most doubtful. As was discussed in the 
previous section, the rotating disks of spiral galaxies show that their baryons carry angular 
momentum and one should expect the dark matter to have similar amounts of angular 
momentum and hence to move on non-radial orbits with distances of closest approach to the 
galactic center at least of order the radius (~ 10 kpc) of the disk. Assumption 4 is motivated 
mainly by simplicity. Below, in the next subsection, we will generalize the model to rid it 
of this assumption. For clarity, we refer to the model with the fourth assumption included 
as the radial infall model. 

Let us call Mj the mass inside at the initial time tj. In a perfectly homogeneous and 
flat universe, Mj is equal to 



Vi = H(ti)fi 



(3.1) 




(3.2) 



if we take tj to be in the matter dominated epoch. Instead, 




2rf 



(3.3) 
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where <5Mj(rj) is a spherically symmetric overdensity. The dark matter shell which is initially 
at radius r, has initially the radial velocity 

Vi{n) = H{ti)n = 2n/3ti , (3.4) 

assuming that £j is small enough so that the very earliest deviations from perfect Hubble 
flow may be neglected. The position r(r,,t) of each shell at time t is determined by solving 
the equations 

d 2 r GM(r,t) 

= — — ' (3 - 5) 



r°° dm- 

M(r,t)= d ri —^e(r-r(n,t)), (3.6) 
Jo dri 

with the initial conditions given in Eq. ( |3.4j) . Q(x) is the Heaviside function, defined earlier. 

The qualitative evolution of the dark matter distribution in phase-space (r, f) may be 
described as follows. Initially, the dark matter particles are located on the line r = H(tj)r. 
As time goes on, this line "winds up" in a clockwise fashion rotating most rapidly near 
r = r = 0. Fig.|I] shows the line on which the dark matter particles are located at a 
particular moment in time. 

The radius r(rj, t) of a given shell initially increases till it reaches a maximum value r*(rj) 
at a time £*(rj). r*(rj) and t*( r i) are called the turnaround radius and turnaround time of 
shell 7*j. After £*(rj), the radius of shell will oscillate with decreasing amplitude. As long 
as it does not cross any other shells, the mass interior to shell is constant, with value Mj, 
and its motion is the well-known motion of a particle attracted by a central mass Mj in the 
limit of zero angular momentum. Shell does not cross any other shells till some time after 
t*(rj), when it is falling onto the galactic center for the first time. Thus, one readily finds: 



r * (r>i)3 and (3.7a) 



2GM l (r l ) 



Mj( ri ) 

After a given shell crosses other shells its motion depends on that of the other shells and 
becomes more difficult to determine. 

Much progress in the analysis of the model came about as a result of the realization 
11,12] that Eqs. (|3.5| ) and (3^) have self-similar solutions for appropriate initial conditions. 



A solution is self-similar if it remains identical to itself after all distances have been rescaled 
by a time-dependent length R(t) and all masses by a time-dependent mass M(t). R(t) is 
taken to be the radius at which dark matter particles are turning around at time t; see 
Fig. |I| M(t) is taken to be the mass interior to R(t) at time t. So, R(t) = r^(rj) and 
M(t) = Mj(rj) with such that t*(rj) = t. A self-similar solution has the properties: 
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M(r,t) = M(t)M(r/R(t)) 



(3.8) 



and 



r{r h t) = r»(ri)\(t/U(ri)) , 



(3.9) 



where M. and A are functions of a single variable. Let us verify that indeed the evolution 
is self-similar for appropriate initial conditions. Substituting Eqs (|3.8|) and (|3~9| ) into Eq. 
fl3.5|) and using Eq. (|3.7a| ), one finds: 



7l~ 



M(t) 



i\ 2 Mi{n) 



M 



(3.10) 



where r = tjt^ivi). We want the RHS of Eq. ( |3. 10| ) to depend only upon r and A(i 
happens for the initial condition: 

SMAn) 



MArt 




This 



(3.11) 



where M and e are parameters, e should be in the range < e < 1, since e = corresponds 
to the extreme case of a r^-independent overdensity whereas e = 1 corresponds to the extreme 
case of an excess point mass located at r = 0. The initial density profile (|3.11|) does not have 
any feature that would distinguish an epoch in the evolution of the galactic halo from other 
epochs. It is this "scale free" property that makes the initial density profile ( |3.11| ) consistent 
with self-similarity, as we are about to show. From now on, for the sake of convenience and 
following Fillmore and Goldreich, we will use Mj instead of T{ to label the shells. Using Eqs. 
( p.3|) and (|3.11|) , and neglecting terms of order 5 Mi /Mi versus terms of order one, we find 
that Eqs. (|3.7|) become 



3tt / AL \ 3e/2 



(3.12a) 



r*(il4 



7T 



2 tl(Mi)GMi 



nl/3 



(3.12b) 



Hence 



At \ 2/3e 
M(t) = M ( — ) , and 
dixt, 



(3.13a) 



R(t) 



8t 2 G 



7l z 



M{t) 



1/3 



(3.13b) 



Therefore 



M(t) 
Mi 



M(t) 



2/3e 



M{U(Mi)) \U(Mi 



T 



2/3e 



(3.14a) 
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2/3+2/9e 



r_JM = R(U(Mj)) = ( U(Mj) ' 
R(t) R(t) \ t 

Thus, Eq. (|3.10|) has the desired form: 

d 2 X 7T 2 t 2 ^ / A 



r -2/3-2/9e 



(3.14b) 



dT 2 8 A 2 \ T 2/Z+2/9e 



M z^Z^7 • (3-15) 



Similarly, using Eq. ( |3.8|) and Eqs. ( |3.14|) , we rewrite Eq. (|3.6| ) as 



M(t) Jo M(t) y w * v " \r*(Mi)JJ 

2 r dr ( A(t) \ 

6 K - ,2/3+2/9. ' ( 3 - 16 ) 



3e 7 1 r l+2/3e 1* r 2/3+2/9e 

which also has the desired form, r varies from 1 to oo. The boundary conditions at r = 1 
are: 

A(l) = l, ^(1) = 0. (3.17) 



Fillmore and Goldreich [|ITJ integrated Eqs. (|3.15 ) and (|3.16| ) numerically for various values 



of e. They also derived analytically the behaviour of -M(£) when ^ — > using adiabatic 
invariants to obtain the motion of the shells. Bertschinger [|12j analyzed the case e = 1. 

Let us mention in passing that the ratio of the density at the turnaround radius to the 
critical density is p*/p c = 97r 2 /[16(3e + 1)]. 



B. Secondary infall with angular momentum. 



As was emphasized earlier, the assumption that the infalling dark matter is devoid of 
angular momentum with respect to the galactic center is inadequate for the calculation of 
the velocity peaks which are the main topic of this paper. However, it is possible to include 
the effect of angular momentum into the secondary infall model while keeping the model 
tractable. We will do this in two steps. First we will assign the same value of angular 
momentum magnitude to all particles in a given shell. Second, we will assign the particles 
in a given shell a distribution of angular momentum magnitudes. The model which results 
from the first step is not so realistic but it is easier to explain. 



1. Single magnitude of angular momentum for all particles on a shell. 

Let's assume first that the particles belonging to shell Mj all have the same magnitude 
of angular momentum Z(Mj) and that they all have, at some initial time, the same radial co- 
ordinate r(Mi,t) and the same radial velocity v r (Mi,t) = dr(Mi,t)/dt. We further assume 
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that at any point f(Mi,t) = fr(Mi,t) on shell Mj, all particles in the shell have their trans- 
verse velocities v±(Mi,r,t) = l(Mi)<p/r(Mi,t) isotropically distributed about f, i. e. each 
direction <p perpendicular to f is equally much represented. As a result of these assumptions, 
each shell remains spherical as it moves through the spherically symmetric mass distribution 
M(r,t) due to all other shells. Initially, at time t,, the shell Mj has radial velocity 

v r (M t ,t) = F(t i )r(M i ,t i ) = 2r( ^ ,tt} . (3.18) 

Eqs. ( |3.2|) and ( |3.3|) hold as before. Provided l(Mi) is not too large, the turnaround radius 
r*(Mj) and time t*(Mj) are still given by Eqs. (|3.7|) to a very good approximation. We will 
use these equations, neglecting the corrections therein due to 1{M{) ^ and find that self- 
similarity is possible. Note, however, that if the corrections to Eqs. ( |3.7| ) due to / (Mj) 7^ are 
included one still finds self-similarity to be possible. The reason we neglect the corrections is 
not to obtain self-similarity but because these corrections are truly small for realistic values 
of angular momentum. 

To obtain self-similarity we assume as a necesary, but no longer sufficient, condition that 
the initial mass distribution is given by the scale free power law of Eq. ( |3 . 1 1[ ) . Eqs. ( |3.12| ), 
( |3.13| ), ( |3.14| ) are then still valid as well. Each dark matter particle satisfies: 



d 2 r _ I 2 GM(r,t) 
dt 2 r 3 r 2 



(3.19) 



Substituting therein 



r (M„ t )=r.(M,)Al__1 , (3.20a) 



M(r, t) = M(t)M U^-J , (3.20b) 

one obtains, using Eqs. ( |3.12|) - ( |3.14| ): 

d 2 \_ l(M l )%(M l ) 2 ( A \ 

dr 2 r,(M,) 4 A 3 8 A 2 ^2/3+2/9, ) . V' ) 

where r = t/t*(Mi) as before. To obtain self-similar solutions, we must make the additional 
assumption that 

W)=i^. (3-22) 

where j is a constant. Then 

<i 2 A ]' 7t 2 T 2 ' 3 >,( X \ ,„ „, 
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Eq. (|3.16|) for -M(£) and the boundary conditions ( p,17| ) remain unchanged. 

The question arises how realistic the model is in the above form. Consider shell Mj near 
its turnaround time t*(Mi). The actual angular momenta of the particles in the shell are of 
course not distributed as in the model. The model assumes v± to be isotropically distributed 
about r. Instead, angular momentum has a unique value l(M{, r) at each point, with Z(Mj, f) 
changing from point to point on the shell. As was discussed in Section |J, r) must have 
at least two zeros on the shell. This implies that there are necessarily some particles in the 
shell which will pass through the center of the galaxy and, by continuity, other particles will 
reach the Earth as well. There are two velocity peaks in the spectrum of cold dark matter 
particles on Earth due to particles falling into and out of the galaxy for the first time, two 
peaks due particles falling into and out of the galaxy for the second time, etc. In contrast, 
in the above model, after a certain galactic age none of the particles falling onto the galaxy 
for the first time reach the Earth because these particles have too much angular momentum. 
Their distance of closest approach to the galactic center exceeds our own distance to the 
galactic center. If we turn for a moment to the results of the computer simulations, we see 
that for the example of Figs.|6| and [7|, in which j = 0.2, only particles which are falling in 
and out of the galaxy for the n-th time with n > 3 can presently reach us. 



2. Distribution of magnitude of angular momenta on a shell. 



In reality, particles at different locations on a given shell have different values of vector 
angular momentum. As a result, the time evolution of a shell is not spherically symmetric 
when angular momentum is present. This is illustrated by Fig.[| in a special axially sym- 
metric case. However, we can restore spherical symmetry by averaging over all possible 
orientations of a physical halo. This corresponds to adopting the model of the previous 
subsection but with a distribution of magnitudes of angular momentum for the particles in 
each shell. Let each shell M, have a fraction rik(Mi) of particles with magnitude of angular 
momentum lk(Mi) where k = 1, . . . , K. To obtain self-similar solutions, n&(Mi) must be 
independent of Mj and 

\2 



rjM, 

The equations for self-similar solutions are then 

a?X k jt it 2 t^.J X 



l ^ =j "im- (3 - 24) 



r k (Mi,t) = n (M l )X k (—. ^— ) , (3.25a) 



dr 2 XI 8 X 2 k \ T 2 / 3 + 2 /^ 



M 37^57 , (3.25b) 



M0 = |X>*f ^e(e-^s) . (3-25C) 
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K 

£ n k = 1 . (3.25d) 

fe=i 

In our calculations, we shall take jk to be distributed according to the density 

^ = !expH7io 2 )- (3-26) 
a J Jo 

Let us explain this choice, starting with the behaviour of dn/dj near j = 0. The angular 
momentum field on a sphere must have at least two zeros. Let us choose the origin (9 = 0) 
of polar coordinates to be at one of them. Assuming that the zero is simple, the Taylor 
expansion of the magnitude of angular momentum function j(9,<p) in powers of 9 starts 
with the term linear in 9: j ~ 9. Then we have: dn/dj ~ dn/d9 ~ 9 ~ j. The cutoff at 
large j in Eq. ( |3.26|) was chosen to be Gaussian for the sake of convenience. The actual 



distribution likely has a sharp cutoff, with a maximum value of angular momentum, but 
that is very similar to a Gaussian. We express our results below in terms of the average over 
the distribution, j = ^Jnj /2. 

A benefit of including angular momentum into the secondary infall model is to produce 
galactic halos with an effective core radius. The radial infall model, i.e. the model without 
angular momentum, produces flat rotation curves for < e < 2/3. Adding angular momen- 
tum has the effect of depleting the inner halo and hence of making the halo contribution to 
the rotation curve go to zero as r — > 0. This is desirable because, in spiral galaxies like our 
own, it is the sum of the contributions from the halo, the disk and the bulge that produces 
flat rotation curves, and the central parts of the galaxy are known to be dominated by the 
bulge and the disk. We define the "effective halo core radius" b as the radius at which half 
of the rotation velocity squared is due to the halo. In our galaxy b is estimated to be of 
order 10 kpc. We will find below that this implies j ~ 0.2 in the model. 

The secondary infall model with a distribution of angular momentum described in this 
subsection still has shortcomings due to the fact that the model averages over all possible 
orientations of a physical halo. In particular, the model is only able to produce estimates of 
the average sizes of velocity peaks. The averages are over all locations at the same distance 
from the galactic center as we are. At some of these locations, a particular velocity peak 
may be much larger than average because that location happens to be close to an inner 
caustic. 



C. Inclusion of baryons 



One may also wish to include the effect of the gravitational potential of the galactic 
bulge and disk. As was already noted, the disk and bulge of our galaxy are conspiring 
with its dark matter halo to produce an everywhere approximately flat rotation curve. We 
may reasonably assume that this was also true in the past because most spiral galaxies are 
observed to have approximately flat rotation curves and they do not all have the same age. 
This suggests a simple way to include the effect of baryons in the self-similar secondary 
infall model, to wit: first obtain for given e the mass function of the model without 
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angular momentum and then use that mass function and Eq. (|3.25b| ) to obtain the phase 
space distribution of the dark matter in the model with angular momentum. We will find 
below that including the gravitational field of the disk and bulge in this manner does not 
have much effect upon the sizes of the highest energy peaks but it does shift their kinetic 
energies upward by deepening the potential well at our location. 



D. e and the spectrum of initial density perturbations. 

The spectrum of the cosmological density perturbations which give rise to galaxies con- 
tains information about the likely value of the parameter e. It has been shown jTJj] that, if 
the density perturbations have a Gaussian distribution, the average density profile around 
a peak in the density distribution is given simply by 

<*(r)>= 5(0)^| (3.27) 

where 5(r) = 5p(r)/p and £(r) =< 5(r)5(0) > is the 2-point correlation function. The latter 
is related to the power spectrum P{k) by 

C(r) = J exp(ik.f)P(k)d 3 k. (3.28) 

The power spectrum is the product P{k) = Ak n T 2 (k) where Ak n is the primordial spectrum, 
taken for simplicity to be a power law, and T(k) is the transfer function. For cold dark 
matter, the transfer function is given by 



= ln(1 + ^' 34g) x [l + 3.89g + (16.1g) 2 + (5.46g) 3 + (6.71g) 4 l 1/4 , (3.29) 
2.34g L J 



where q = kMpc/h 2 . The Harrison-Zel'dovich spectrum corresponds to n = 1. Eqs.( p.27|) 



and ( |3.28| ) imply that if P{k) ~ k a on some momentum scale k, then £(r) ~ r - Q - 3 anc l hence 
e = (a + 3)/3 on the corresponding length scale r = 1/k. We computed a = dlnP/dlnk for 
n = 1 and plotted the resulting e(r) in Fig. [4] for the relevant scales. The figure suggests 
that e is of order 0.2 - 0.3 on the galactic scale. 



IV. NUMERICAL INTEGRATION 



In this section, we present the results from numerically integrating Eqs. (|3.25b| , p!25c|) 
for various values of the parameter e and various angular momentum distributions, includ- 
ing no angular momentum, a single value of angular momentum and the distribution of 
Eq.( |3.26| ). The function A(r) gives us the phase-space distribution of the particles through 
the equations: 



r{Mi,t) =r,(Mi)X 



R(t)X(r)T 



-2/3-2/9e 



(4.1a) 
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v{M h t) = = M r l/3-2/9^. (41b) 

at t dr 

If there is a distribution of angular momentum values, the functions A(r), r(M i) t) and 
v(Mi,t) carry an index k which we have suppressed here to avoid cluttering the equations. 

To solve Eq. fl3.25rj|) for the particle trajectory A(r) we need to know the mass function 



-M(£). This function, in turn, is given in terms of the trajectory A(r) by Eq. (|3.25c|) or, 
equivalently, by: 

M(0 = £ (r^CO - r 2 - n 2/3t (0) , (4-2) 

n=l 

where the Tj(£) correspond to the moments of time when the trajectory crosses radius 
r = £R(t), i.e. they are the solutions of A(r) = ^ r 2 / 3 + 2 / 9e . Following Fillmore and Goldreich 
|Ti|| , we solve Eqs.( |3.25b~D and ( |3.25cj ) simultaneously by a technique of successive iterations. 
Starting with some arbitrary mass profile -M(£) (we took -M(£) = £ 2 ) we find A(r), which 
is then used to derive a new mass profile, from which a new trajectory is derived, and so 
on. The procedure is repeated till it converges. We find that the mass profile changes very 
little after 5 iterations. Typically we run 10 iterations to get the final results. 

Fig. [5] shows the phase-space diagram for the case e = 0.2 and zero angular momentum. 
The solid line in that figure shows the location of all the particles in phase space at a 
given time, i.e. it is the set of points (r(M i ,t),v(M i ,t)) for all Mj. The radial distances 
are normalized to the turnaround radius R at time t and the velocities are normalized to 



GM(t)/R(t) = rcR(t)/\/8t which is the rotation velocity at the turnaround radius. Fig. || 
shows the phase space diagram for the case e = 0.2 and a single value of angular momentum 
j = 0.2. The particle trajectory A(r) for that case is shown in Fig. |7|. 

A convenient way to show the mass distribution is by showing the rotation curve. We 
define z/(e, £) by: 

iH^MM/r^j^jJ. (4.3) 

With this definition, we have z/(e,£ =1) =1- The functions u(e, £) 2 obtained by numerical 
integration are plotted in Fig. |8| for various values of e and j = 0. 

To fit the model to our galactic halo, we must choose values of the present turnaround 
radius R = R(t) and of M = M(t). Equivalently, we may choose values of the present age t 
and of R. M is given in terms of R and t by Eq. (|3.13b| ). t is given in terms of the Hubble 



rate H = hlOO km s _1 Mpc -1 by the relation t^ 1 = 3H /2. We will use h to state the 
age of the universe. Then we fix R in terms of h by requiring that the model reproduce the 
measured value, v rot = 220 km s _1 , of the rotation velocity in our galaxy. Let us call u(e) 
the value of i>(e, £) in the flat part of the rotation curve, near r = 0.02R for e < 0.4; see Fig. 
|l z/(e) is related to v rot by Eq. ([4.3|). This implies: 



Rh = 1.32z/(e)- 1 Mpc. (4.4) 
Table I gives z^(e) 2 , Rh and Mh for various values of e. 



18 



Our distance to the galactic center is taken to be 8.5 kpc and we define £ s = 8.5 kpc/R. 
The contribution of the n-th velocity peak to the local halo density is given by: 



Pn 4K&a | 9e r n A n -(66 + 2)A n |' ( ' j 

where A = dX/dr, and A ra and r n are the solutions of A(r) = ^ s r 2 ^ 3+2 ^ 9e . The kinetic energy 
(in a frame of reference which is not rotating along with the galactic disk) per unit particle 
mass in the n-th peak is given by: 

E = -v 2 = — T 2/3 - i/9€ (— + \ 2 \ (4 6) 

n 2 " 2t 2 n \K ) ' 1 ' 

We shall express E n in units of (300kms~ 1 ) 2 /2 when presenting our results. 

We now discuss in greater detail the results specific to the different types of angular 
momentum distributions used. 



A. Radial infall 



Without angular momentum all particles pass through the origin, r = 0, at each oscilla- 
tion. To avoid infinities in the numerical integration, a regulator at small r is required. The 
one which is most convenient for us and which we use is to give a small amount of angular 
momentum to the dark matter particles. We found j 2 = 10~ 6 to be small enough for our 
purposes. 

Fig. |5] shows a typical phase space distribution. Fig. [8] shows the rotation curves 
for various values of e. An analytical treatment of the radial infall model using adiabatic 
invariants predicted Jill the behaviour of the density near the origin to be: p oc r - 9e /( 3e + 1 ) 
in the range 2/3 < e < 1 and p oc r -2 in the range < e < 2/3. These predictions agree 
very well with our results. The rotation curves do indeed go to a constant near the origin 
when < e < 2/3 except for small logarithmic corrections. The analytical treatment given 
in section V suggests the behaviour p ~ l/(r 2 -y/ln(l/r)) for small e. 

Fig. shows the velocity peaks on Earth predicted by the radial infall model with e = 0.2 
and h = 0.7. The rows labeled j = 0.0 in Table |lj give the density fractions and kinetic 
energies of the five most energetic incoming peaks for the cases e = 0.2 and 1.0, and h = 0.7. 
For each incoming peak there is an outgoing peak with approximately the same energy and 
density fraction. We find that, in the radial infall model, the sizes of the two peaks due to 
particles falling in and out of the galaxy for the first time are large, each containing of order 
10% of the local halo density for e in the standard CDM model inspired range of 0.15 - 0.4. 
As was emphasized already, this spectrum is unrealistic because angular momentum has a 
non-negligible effect upon the peak sizes. 
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B. Single non-zero value of angular momentum. 



Fig. ^ shows the phase space diagram for the case e = 0.2 and j = 0.2. In the model 
with a single non-zero value of angular momentum, the halo distribution has a set of inner 
caustics in addition to the usual outer caustics. However, the inner caustics are spheres in 
the model whereas they are rings for a physical halo. The density profile for the model is 



shown in Fig. [T(J. There is a break in the logarithmic slope near the first inner caustic, near 
r = 0.01.R in the figure. For e < 2/3 we find p oc r~ 2 outside the first inner caustic (the same 
as with j = 0) but p oc r~ 7 inside with 7 = 9e/(3e + 1). This observation is related to the 
fact that we find the function A(r) to be oscillating with constant amplitude and constant 
period for large r when j 7^ 0; see Fig. 0. This can happen only if t 2 ^ 3e M. (A/t 2 / 3+2 / 9<: ) in 
Eq. ( |3.25b| ) is independent of r at large r. For M.{r) oc r a , this implies a = 3/(3e + 1) or 
7 = 3 — a = 9e/(3e + 1). We find that the exponent a does not depend upon j, but the 
radius r c at which the break in power law behaviour occurs does depend upon j: the smaller 
j, the smaller r c . 

These observations may be understood as follows. When angular momentum is present, 
the contribution of a single phase-space sheet to the density profile p(r) is non-singular near 
r = 0. In that result of self-similarity the mass profile M(r), at r small enough 

that the influence of the first one or two phase-space sheets may be neglected, has the same 
functional dependence upon r as the dependence of M upon R after eliminating t from the 
expressions for M{t) and R(t) in Eqs. (|3.13|) . This yields M oc r 3 ^ 3t+1 \ The fact that p(r) 
behaves as a negative power of r near r = is in agreement with the results of N-body 
simulations [19]. The e-dependence of this power law was already known to White and 
Zaritsky |IJ. 

The spectrum of velocity peaks that the model with a single value of angular momentum 
predicts is unrealistic. In particular, there are in this model no peaks on Earth associated 
with particles falling in or out of the galaxy for the n-th time with n small, because such 
particles have too much angular momentum to reach the Earth radius. In the case of Fig. 
H there are only peaks for n > 3. But as we argued at length in section II, there are in 
reality peaks on Earth due to particles falling in and out of the galaxy for the n-th time 
with n = 1, 2, 3, ... Angular momentum reduces the sizes of the peaks with small n but does 
not suppress them completely. 



C. Distribution of angular momenta. 



This model describes a physical halo averaged over all possible orientations. The particles 
are assumed to have the distribution of angular momenta given in Eq.( [3.26| ). For the purpose 
of numerical integration, the angular momentum was discretized with a spectrum of four 
hundred values. The phase space diagram is like the one of Fig. [| except that there are four 
hundred such diagrams superimposed on one another, one for each value of j. 

As was mentioned earlier, angular momentum has the effect of making the contribution 



of the halo to the rotational velocity go to zero at the galactic center as shown by Fig. Ill 



We define the "effective core radius" b to be such that z/ 2 (e, b/R) = u 2 (e)/2. We find that 
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near r = 0, the density p(r) oc r~ 7 with 7 = 9e/(3e + 1) as in the case of a single value 
of angular momentum. This behaviour is expected for the same reason as we gave for that 
case. The only change is that the transition between the region at large r where p{r) oc r~ 2 
and the region at small r where p(r) oc r~ 7 is smoother now because the critical radius r c 
where the transition occurs depends upon j and there is now a distribution of j values. 

The velocity peaks for j = 0.2, e = 0.2 and h = 0.7 are shown in Fig. |1^. The velocity 



peaks for the same choice of parameters except j = 0.4 are shown in Fig. |13| to indicate the 
sensitivity of the peaks upon j. The spectrum of velocity peaks is also sensitive to the value 
of e. It is shown for the cases e = 0.15 and e = 1 on Figs. [14| and [15| respectively. 

Table |I| gives the values of the current turn-around radius R, the effective core radius b, 
the halo density at our location p(r s ), and the density fractions and kinetic energies of the 
five most energetic incoming peaks for various values of e, j and h. 



D. Self-similar infall with baryons 



As was discussed earlier, the dark matter phase space distribution for the input param- 
eters e, h and j is obtained in this case by solving Eq. ( 3.25cj ) using the mass distribution 
-M(£) for the same values of e and h but j = 0. The resulting velocity peaks are shown in 
Fig. [16] for e = 0.2, h = 0.7 and different values of j. Since the particles for different values 
of j but the same values of e and h are all moving in the same gravitational potential, the 
kinetic energies per unit particle mass E n of the peaks depend only very weakly upon j. 
So, we have combined the spectra for different j values in one figure. Note that the E n are 
larger in this case than in the case of Figs. [12] and because the gravitational well near 
the center of the galaxy is deeper. Fig. O shows the results for e = 0.4, h — 0.7. 



V. AN ANALYTICAL TREATMENT. 



In this section, we derive analytical expressions for the sizes and locations of the velocity 
peaks. The treatment involves the following approximations: 

1. the mass distribution, including the contributions from both baryons and dark matter, 
is taken to be M(g) = £ for < £ < 1 

2. the dimensionless angular momentum j values are assumed to be small 

3. our distance r s to the galactic center is assumed to be small compared to the oscillation 
amplitudes of the particles in the peaks under consideration. 

The method of adiabatic invariants will be used to obtain the motion of the dark matter 
particles. 

We first treat the case of a single angular momentum value, i.e the model described in 
subsection III.B.l. Since (£) = £, the equation of motion in rescaled variables of a particle 
with angular momentum j is: 
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d 2 X j 2 IT 2 2,2 
— r 3^ _1 J 



dr 2 A 3 



IX 



(5.1) 



Let \m and X m be respectively the amplitude of oscillation and the distance of closest 
approach at time r. In the spirit of the method of adiabatic invariants, Am and X m are 
assumed to be slowly varying functions of r. This is a valid assumption for all oscillations 
except the first one. The first oscillation is only marginally adiabatic. For given Am, the 
velocity at time r and position A is: 



dX _ _ 
dr 
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The adiabatic invariant is: 
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In the limit of small j, y 22 - may be neglected and we obtain Xm(t) 2 t 
Since Am(1) = 1, we have: 



constant. 



7- ^ 



(5.4) 



The period of oscillation is: 
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Using Eq. ( |5.4| ), we have in the limit of small j: 
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Let us introduce the "phase" <f(r) : 
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dr 
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Using Eq. (|5.6|) and setting y?(l) = 0, we have: 

(p{r) = 9e 



7T V ~ T3(3t +1 ) — 1 



4 4 + 3e 

The times r n at which the particle passes by the solar radius r s are given by: 



(5.8) 
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7T 



<p( Tn ) = -(2n-l) + o(r./R) 



(5.9) 



There are two peaks for given n, one ingoing and the other outgoing. The differences between 
the properties of the two peaks are of order o(r s /R) and are neglected. 

The contribution of each of the two n-th velocity peaks to the local density is given by: 



Pn 



M 



6r n 3 



2n — 2.1 

3\ x 3e> 



A-nRrl |9er n f ( r „) - (6e + 2)A(r n 



(5.10) 



where 



K T n) = Tn 



2,2 
3^96 



(5.11) 



and jr(T n ) is given by Eqs. Q5.2|) and ( |5.4j ). We neglect the second term in the denominator 
of Eq. ( |5.10| ) since it is o(r s /R) relative to the first term. Combining everything, we have 
the following estimate for the peak sizes due to dark matter particles with a single value of 
angular momentum j: 
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where 
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(5.13) 



Note that p n is the local energy density contributed by each of the two n-th peaks. When the 
expression under the square root in Eq. (|5 . 1 2 ) is negative, one must set the corresponding 
p n = since this corresponds to the situation where the particles have too much angular 
momentum, and hence too large a distance of closest approach to the galactic center, to reach 
the solar radius. The kinetic energy per unit particle mass of the dark matter particles in 
the n-th peaks is: 



E,, 



J 2 
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R 






JsTn 



+ o{r s /Rf 



(5.14) 



R and t are determined in terms of h by t = 2/3Hq and Eq. ( |4.4| ) as before. 

For the case e = 0.2, j = 0.2, h = 0.7, the quantity under the square root in Eq. ( |5.12| ) 
is negatve and hence p n = for n=l,2 and 3. This is consistent with the phase space 
distribution shown in Fig. |B] which was obtained by numerical integration. We found Eqs. 
( |5.12j - I5.14Q to be consistent with the results from numerical integration at the 30% level 
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or so. The agreement is worse for the case of zero angular momentum, probably because 
the logarithmic singularity of the potential at the origin in that case makes the motion non- 
adiabatic there. Also, when comparing results for the E n values, one should allow for an 
overall shift between the two calculations due to a change in the depth of the gravitational 
potential at the solar radius. Indeed, the analytical calculation has the mass distribution 
= £ which implies a perfectly flat rotation curve, whereas the rotation curves for the 
numerical calculations £1X6 cLS shown in Fie.S. 



The formalism readily accomodates a distribution dn/dj of angular momenta. Eq. (|5.14|) 
for the peak kinetic energies still applies. The expected peak sizes, in the sense of an average 
over all locations at distance r s from the galactic center, are given by the convolution of 
dn/dj with the expression, Eq. ( |5.12j ), for the peak sizes when there is a single value of 
angular momentum. Thus: 
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where r n is given by Eq. ( 5.13| ) as before and j n is the maximum value of j for which the 
argument of the square root is positive: 
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For the angular momentum distribution of Eq. (3.26), one has 
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where 



F(u) 
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dxe~ 
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A graph of F(u) is shown in Fig. [18 . Table 3 shows the peak sizes predicted by Eq. flS.lTp for 
the first eight peaks in the case e = 0.2, j = 0.2 and h = 0.7. The peak sizes agree with the 
results of the numerical integration, given in the fifth line of Table 3, to within 15%. The 
peak energies also agree within 15% after one has subtracted an overall shift caused by the 
fact that the gravitational potential at our location is considerably deeper in the analytical 
calculation than in the numerical one. 



VI. CONCLUSIONS 



Motivated by the prospect that the spectrum of dark matter particles on Earth may 
some day be measured in a direct detection experiment, we have endeavoured to obtain 
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predictions for the properties of that spectrum. Previously, it had generally been assumed 
that the spectrum is isothermal. In contrast, we find that there will be large deviations from 
an isothermal spectrum in the form of peaks in velocity space associated with particles that 
are falling in and out of the Galaxy for the first time and with particles that have fallen in 
and out of the Galaxy only a small number of times in the past. 

To obtain estimates for the velocity magnitudes of the particles in the peaks and for 
the contributions of the individual peaks to the dark matter local density, we have used the 
secondary infall model of galactic halo formation. We have generalized the extant version 
of that model to include the effect of angular momentum. We forced spherical symmetry 
onto the model with angular momentum by averaging over all possible orientations of a 
physical halo. As a result, the model can only make predictions for the average properties 
of the velocity peaks, the average being over all locations at the same distance from the 
galactic center as we are. We found that the model with angular momentum has self-similar 
solutions if the angular momentum distribution, as well as the initial overdensity profile, has 
a particular scale-free form. The self-similar solutions were analyzed in detail numerically 
and analytically. 

The model produces a good overall fit to what is known about galaxies like our own. It 
produces flat rotation curves when the parameter e is in the range to 2/3. The expected 
value of that parameter in models of large scale structure formation with cold dark matter 
and a flat (Harrison- Zel'dovich) spectrum of primordial density perturbations is e ~ 0.25. 
The model implies a relationship (cfr. Eq. (|4.4j ) and Table I) between the current turn- 
around radius R, the present age of the universe, the galactic rotation velocity and the 
paramemter e. This relationship is consistent with observations and e ~ 0.25. 

In the model the galactic rotation curve is approximately flat all the way out to the 
turn-around radius R. R is of order 1-2 Mpc for our galaxy. So far, the rotation curves of 
individual galaxies have been measured, and have been found to be flat, up to distances of 
order 100 kpc, implying masses of order 10 12 M Q or larger. The discovery of flat rotation 
curves |21[] caused a well-known revolution in our concept of galaxies. Prior to these mea- 
surements galaxies were thought to have size of order 10 kpc and mass of order 10 ll M Q . 
The model implies that galaxies like our own are of order 20 times bigger even than the 
minimum size implied by the rotation curve measurements. The size is 1.8 Mpc and the 
mass is 1.7 x 1O 13 M for the case e = 0.25, h = 0.7. See Table I for other cases. With 
galaxies that massive, one has Q ~ 1 in galaxies. 

We found that the effect of angular momentum is to deplete the inner parts of the halo 
with the result that the halo contribution to the rotation curve goes to zero at r = 0. 
The model establishes a direct relation between the amount of angular momentum and the 
effective core radius b, defined as the radius at which the halo contributes half of the rotation 
velocity squared. Table [H] gives the turn-around radius R, the effective core radius b and 
the local halo density p as a function of the input parameters: e, the average amount of 
dimensionless angular momentum j and the Hubble parameter h. The observed value (220 
km/s) of our galaxy's rotation velocity and our distance to the galactic center (8.5 kpc) were 
used as fixed input parameters. 

Table [H] also gives the average values of the peak sizes as fractions of the local halo 
density p and the kinetic energies per unit particle mass of the particles for the first five 
incoming peaks as a function of the variable input parameters e, j and h. For each incoming 
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peak there is an outgoing peak with approximately the same kinetic energy and average local 
density. Let us emphasize again that the peak sizes given are averages over all locations at 
the same distance (8.5 kpc) from the galactic center as we are. It is not possible to make 
more precise predictions for the peak sizes on Earth without assuming a particular angular 
momentum field for the infalling dark matter. 
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TABLES 
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TABLE II. Density fractions f n and kinetic energies per unit particle mass E n of the first five 
incoming peaks for various values of e, h and the average dimensionless angular momentum j. Also 
shown are the current turn-around radius R in units of Mpc, the effective core radius b in kpc, 
and the local density p in units of 10 -25 g cm -3 . The /„ are in percent and the E n are in units of 
0.5 x (300 km s" 1 ) 2 . 
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TABLE III. 


Values of j n , E n 
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FIG. 1. Phase space distribution of the halo dark matter particles at a fixed moment of time. 
The solid lines represent occupied phase-space cells. The dotted line corresponds to the observer 
position. Each intersection of the solid and dotted lines produces a velocity peak. 




FIG. 2. Positions in physical space at successive moments in time t\ < t<i < ■ ■ ■ < t% oi the 
particles on a turnaround sphere that is intially rotating rigidly about the vertical axis. The * 
indicates the appearance of an inner caustic ring. 
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log(r) 

FIG. 3. A small scale sub-clump falling into the galaxy for the first time. 
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FIG. 4. The expected value of the e parameter as a function of galaxy size, in models of structure 
formation based upon cold dark matter and a flat (Harrison-Zel'dovich) spectrum of primordial 
density perturbations. We defined /jq.7 = h/0.7 
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FIG. 5. The phase space distribution of halo dark matter particles at a fixed moment of time 
for the case e = 0.2 and j = 0. The solid lines represent occupied phase space cells. The dotted 
line corresponds to the Sun's position if h = 0.7. 
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FIG. 7. The function A(r) for e = 0.2, j = 0.2. 
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FIG. 8. Rotational velocity squared curves for different values of e and j = 
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FIG. 9. The density fractions f n = p n /p (in percent) and the kinetic energies per unit particle 
mass E n of the peaks at the Sun's position, for e = 0.2, j = and h = 0.7. The peaks form pairs. 
One member of each pair is due to particles with positive radial velocity and the other is due to 
particles with negative radial velocity 
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FIG. 10. Density profile for the case e = 0.2, h = 0.7 and a single value of angular momentum 
.7 = 0.2. 
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FIG. 11. Rotational curves for the case e = 0.2, with and without angular momentum. 
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FIG. 12. The same as Fig. | but e = 0.2, h = 0.7 and j = 0.2. 
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FIG. 13. The same as Fig. | but e = 0.2, h = 0.7 and j = 0.4. 
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FIG. 14. The same as Fig. | but for e = 0.15, h = 0.7 and j = 0.2. 
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FIG. 15. The same as Fig. | but e = 1, h = 0.7 and j = 0.2. 
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FIG. 16. The same as Fig. ^ but including the contribution of baryons to the galactic gravita- 
tional potential, e = 0.2 and h = 0.7 in all cases. The peaks are shown explicitly for j = 0.1. The 
dashed lines go through the tops of the peaks for the cases j = 0.2 and j = 0.4 while the dotted 
line corresponds to j = 0. 
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